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The  parabolic  approximation  method  is  widely  recogniz 
and  predicting  sound  transmission  intensity  in  diverse 
attractiveness  is  that  solutions  are  marched  in  range,  thi 
storage  requirements  when  using  the  full  wave  equatioi 
implementations  employ  a  range  step  size  that  is  presc 
that  remains  fixed  for  the  duration  of  the  computation 


as  useful  for  accurately  analyzing 
n  environments.  One  reason  for  its 
eby  avoiding  the  large  internal 
Present  finite-difference 
bed  by  either  the  user  or  the  code  and 
n  algorithm  is  presented  in  which  the 


range  step  is  adaptively  selected  by  a  procedure  within  a  version  of  the  implicit  finite-difference 
( IFD )  implementation  of  the  parabolic  approximation.  An  error  indicator  is  computed  at  each 
range  step,  and  its  value  is  compared  to  an  error  tolerance  window  that  is  readily  specified  by 
the  user.  If  the  error  indicator  falls  outside  this  window,  a  new  range  step  size  is  computed  and 
used  until  the  error  indicator  again  leaves  the  tolerance  window.  Furthermore,  for  a  given 
tolerance,  this  algorithm  generates  a  range  step  size  that  is  optimal  in  a  specified  sense  and  that 
often  leads  to  large  decreases  in  run  time.  Additional  related  modifications  to  the  IFD 
implementations  are  discussed.  Several  examples  are  presented  that  illustrate  the  efficacy  of  the 


enhanced  algorithm,  r, 
PACS  numbers:  43  30.Bp 
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INTRODUCTION 

Since  its  introduction  to  the  underwater  acoustics  com¬ 
munity  more  than  a  decade  ago,'  the  parabolic  approxima¬ 
tion  to  the  reduced  wave  equation  has  seen  widespread  appli¬ 
cation  in  many  important  underwater  sound  propagation 
problems.  This  type  of  approximation  easily  and  accurately 
handles  strongly  range-dependent  environments,  offering  a 
capability  beyond  the  reach  of  normal  mode  theory  and  re¬ 
lated  methods.  Since  it  is  effective  at  low  frequencies,  it  is 
suitable  for  analyzing  and  predicting  propagation  in  situa¬ 
tions  where  diffraction  may  be  important,  thus  offering  a 
significant  advantage  over  ray  theory.  Another  reason  for 
interest  in  this  method  is  computational  efficiency.  The 
fruits  of  this  approximation  are  parabolic  equations  (PEs) 
that  can  be  marched  in  range,  thereby  saving  considerable 
amounts  of  internal  core  storage. 

Numerical  implementations  of  PEs  have  been  devel¬ 
oped  using  different  algorithms.  One  of  these,  using  an  im¬ 
plicit  finite-difference  scheme,2  is  well  suited  to  ocean  envi¬ 
ronments  in  which  there  are  significant  amounts  of  acoustic 
interaction  with  the  ocean  bottom.  Other  algorithms,  based 
on  the  split-step  method3  (which  uses  fast  Fourier  trans¬ 
forms),  have  proven  advantageous  for  modeling  deep-ocean 
sound  channel  propagation.  Both  types  of  algorithms  have 
been  refined  in  many  ways.  For  example,  the  implicit  finite- 
difference  algorithm  has  been  extended  to  handle  wide-angle 
PEs,4  while  the  split-step  method  has  been  implemented  on  a 
high-speed,  special-purpose  array  processor.3 

Although  PEs  offer  significant  computational  advan¬ 
tages  over  methods  for  the  (elliptic)  Helmholtz  equation, 
many  propagation  problems  pose  formidable  computational 


problems  and  may  be  infeasible  with  present  codes.  For  ex¬ 
ample,  results  in  the  time  domain  can  be  obtained  by  solving 
a  PE  for  many  different  frequencies  in  the  same  channel  and 
then  superposing  the  results  with  Fourier  transforms.6 
Broadband  problems  solved  by  this  method  require  dozens, 
perhaps  hundreds,  of  computations  at  different  frequencies. 
Also,  high-frequency  propagation  problems  may  require 
computational  meshes  of  prohibitively  small  size,  making 
solutions  simply  too  expensive  to  compute.  Deep-ocean 
channels  of  long  range  can  require  large  amounts  of  time, 
limiting  the  ability  to  perform  simulated  propagation  studies 
under  a  variety  of  interesting  conditions.  In  addition,  certain 
classes  of  three-dimensional  problems  can  be  computed  by 
two-dimensional  solutions  through  vertical  slices  of  the 
ocean.7  For  sufficiently  large  azimuthal  regions,  this  method 
clearly  requires  large  amounts  of  time  for  each  computation. 

There  is  a  clear  and  pressing  need  for  substantially  faster 
propagation  codes."  In  part,  this  requirement  for  speed  will 
be  satisfied  with  more  powerful  computers,  such  as  array 
processors  and  supercomputers,9  which  may  become  avail¬ 
able  in  a  few  years.  Even  so,  computational  acoustic  algo¬ 
rithms  must  be  made  as  efficient  as  possible,  regardless  of  the 
type  of  machine.  One  way  to  significantly  contribute  toward 
the  goal  of  obtaining  optimal  efficiency  is  through  the  use  of 
adaptive  computational  methods.  In  these,  the  mesh  consist¬ 
ing  of  the  set  of  discrete  points  at  which  the  solution  is  com¬ 
puted  is  adjusted  so  that  the  number  of  points  used  is  mini¬ 
mized  while  some  measure  of  the  solution  accuracy  is  pre¬ 
served.  For  example,  one  implementation  of  the  split-step 
method*  employs  adjustments  or  range  step  size  based  on 
truncation  error  estimates  that,  in  turn,  depend  on  partial 
derivatives  of  the  refractive  index.  In  addition,  depth  step 
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size  can  increase,  but  not  decrease,  in  response  to  specialized 
conditions  obtained  from  spectral  energy  estimates. 

It  is  the  purpose  of  this  paper  to  describe  a  method  for 
accomplishing  part  of  this  objective  for  finite-difference 
methods.  We  emphasize  here  that  the  approach  we  describe 
is  substantially  different  from  the  one  employed  in  Ref.  3.  In 
addition,  the  focus  of  our  work  is  on  the  performance  accu¬ 
racy  on  the  adaptive  algorithm  as  well  as  its  efficiency.  In 
Sec.  I,  a  two-dimensional  narrow-angle  parabolic  approxi¬ 
mation  and  a  finite-difference  algorithm  implementation  for 
solving  the  resulting  PEs  are  reviewed.  Notation  pertinent  to 
the  discretized  PE  solution,  as  well  as  a  discussion  of  current 
methods  of  selecting  fixed  mesh  sizes,  are  presented.  An  er¬ 
ror  indicator  is  introduced  and  a  method  for  adaptively  se¬ 
lecting  the  range  step  using  this  indicator  is  discussed  in  Sec. 
II.  In  Sec.  Ill  are  presented  several  example?  that  demon¬ 
strate  the  computational  advantages  attained  by  adaptively 
selecting  the  range  step.  Finally,  Sec.  IV  contains  a  summary 
of  our  principal  results. 

I.THE  PARABOLIC  APPROXIMATION  AND  IFD 
IMPLEMENTATION 

Parabolic  approximations  typically  begin  from  the 
Helmholtz  equation 

V2p  4-  kln2p  =  0,  (1) 

which  governs  the  acoustic  pressure  field  p  in  a  steady,  quies¬ 
cent  medium  due  to  a  harmonic  point  source  of  frequency  / 
In  Eq.  ( 1 ),  k0  =  2irf  /c„  is  the  wavenumber,  n  =  c(/c(r,<?,z) 
is  the  index  of  refraction,  c(  r,6j)  is  the  sound  speed,  and  c()  is 
a  reference  sound  speed.  The  farfield  assumption  is  given  by 

p  =  u(r,e^)H{0'\k0r),  (2) 

where  H(0n(k0r)  is  the  Hankel  function  of  the  first  kind  of 
order  zero.  Equation  (2),  assumed  to  hold  for  sufficiently 
large  values  of  r,  supposes  that  there  are  only  outward  travel¬ 
ing  waves,  i.e.,  no  backscattering.  The  quantity  u(r,Qj)  is  a 
slowly  varying  function  of  position  that  modulates  the  Han¬ 
kel  function.  If  Eq.  (2)  is  substituted  into  Eq.  ( 1 ),  and  ap¬ 
propriate  conditions  are  satisfied,  the  following  “standard 
PE”  can  be  shown  to  hold: 

lik0ur  +  u„  +  k  l(n2  -  1)«  =  0.  (3) 

Since  no  6  derivatives  appear  in  Eq.  (3),  we  suppress  the  6 
dependence  of  u  and  write  u(rj).  Details  on  the  derivation 
of  Eq.  (3)  can  be  found  in  Refs.  I,  10,  and  elsewhere.  It  is 
important  to  note  that  Eq.  (3)  is  not  the  only  parabolic  ap¬ 
proximation  that  could  be  derived.  Other  examples  include 
PEs  for  three-dimensional  propagation"  and  wide-angle 
propagation,12  and  additional  ones  can  be  obtained  from 
more  general  versions  of  Eq.  ( 1 ) .  For  example,  PEs  have 
been  derived  that  are  appropriate  for  sound  channels  in 
which  the  medium  is  moving1"  or  the  density  is  variable. 

There  are  several  numerical  algorithms  available  for 
solving  Eq.  (3),  together  with  an  appropriate  set  of  bound¬ 
ary  and  initial  conditions.  In  particular,  one  using  implicit 
finite  differencing  of  the  partial  derivatives  has  been  widely 
distributed  in  the  underwater  acoustics  community.  Details 
of  this  finite-difference  algorithm  can  be  found  in  Ref.  13,  It 
utilizes  a  Crank-Nicolson  scheme  to  march  the  solution  in 


range  with  a  formal  order  of  accuracy  of  Of  (Ar) 2  +  Afz2)], 
where  Ar  and  A z  are  the  range  and  depth  step  sizes,  respec¬ 
tively.  This  method  is  popular  for  parabolic  systems  because 
of  its  accuracy  and  absolute  stability.  One  advantage  to  this 
approach  is  its  particularly  convenient  ability  to  model  hori¬ 
zontal  as  well  as  irregular  interfaces  between  layers  with  dif¬ 
ferent  acoustic  properties.  The  implementation,  known  as 
implicit  finite-difference  (IFD)  implementation,  contains  a 
number  of  powerful  and  flexible  features  that  permit  its  ap¬ 
plication  to  a  broad  variety  of  ocean  acoustic  environments. 
For  instance,  an  input  stream  provides  a  mechanism  for  the 
user  to  prescribe  source  data,  specify  selected  numerical  pa¬ 
rameters  such  as  mesh  size,  describe  many  types  of  complex 
propagation  environments,  and  select  a  variety  of  output  op¬ 
tions.  IFD  implements  pressure-release  and  rigid  boundary 
conditions  on  flat  or  sloping  boundaries,  and  matches  the 
wave  field  from  the  water  column  into  any  number  of  sub¬ 
bottom  layers  with  different  densities,  sound  speeds,  and  vol¬ 
ume  attenuations.  It  also  contains  a  feature  to  apply  an  artifi¬ 
cial  absorbing  layer  (beneath  sub-bottom  layers)  that  is  used 
to  enforce  a  pressure-release  bottom  boundary  in  many  ap¬ 
plications.  The  code  generates  several  output  files,  including 
ones  containing  the  complex-valued  solution  to  the  PE,  as 
well  as  the  transmission  loss.  Additional  features  are  avail¬ 
able  and  are  documented  in  Refs.  2  and  13.  For  the  rest  of 
this  paper,  IFD  will  refer  to  the  implicit  finite-difference 
implementation  described  in  Ref.  2  (and  modified  by  us  for 
increased  speed  and  accuracy),  while  EIFD  will  refer  to  the 
implementation  including  the  enhancements  discussed  here¬ 
in. 

In  the  IFD  implementation,  the  numerical  solution  of 
Eq.  (3)  is  represented  as  u ",  where 

u”  =  u(rm,z„),  n  =  1,...,.V.  (4) 

In  Eq.  (4),  (rm,z„ )  corresponds  to  a  point  on  the  range- 
depth  mesh,  and  M  and  N  are  integers  indicating  the  maxi¬ 
mum  number  of  range  and  depth  points  on  the  mesh.  Al¬ 
though  the  algorithm  treats  nonhorizontal  boundaries  by 
adding  or  deleting  depth  mesh  points,  we  take  N  fixed 
throughout  this  paper.  We  note  that  this  version  of  IFD  only 
handles  horizontal  interfaces,  so  that  sloping  or  irregular 
interfaces  are  approximated  by  stepwise  horizontal  inter¬ 
faces. 

IFD  requires  users  to  either  specify  the  characteristic 
dimensions  of  the  computational  grid  upon  which  the  solu¬ 
tion  is  computed  or  to  accept  default  values.  With  IFD,  A r 
and  Az  represent  fixed  range  and  depth  increments,  so  that 
rm  =  mAr  and  z„  —  n\z.  It  is  not  always  clear  a  priori 
whether  there  are  optimal  choices  for  A r  and  Az  dimen¬ 
sioned  to  minimize  computation  time  while  preserving  ap¬ 
propriate  measures  of  solution  accuracy.  In  the  standard 
IFD  implementation,  the  typical  default  values  are  Ar  =  A  / 
2  and  Az  =  A  /4,  where  A  is  the  acoustic  wavelength  of  the 
source  signal.  In  general,  these  mesh  sizes  tend  to  be  un¬ 
necessarily  small.  Normally,  this  can  only  be  corrected  when 
the  user  has  had  experience  with  computational  results  for  a 
particular  problem. 

The  geometry  of  the  sound  channel,  together  with  val¬ 
ues  of  all  important  parameters,  is  provided  to  IFD  via  an 
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input  runstream  file.  Our  implementation  adds  two  addi¬ 
tional  parameters  to  this  runstream:  an  error  tolerance  e  and 
a  “cut-in”  range  value.  These  parameters  will  be  discussed  in 
detail  in  Sec.  II.  For  output,  IFD  generates  files  containing 
the  complex  valued  solution  u(rm,zn )  and  the  transmission 
loss.  Our  implementation  appends  three  additional  files  that 
contain  the  step  size  used  at  each  range,  an  indicator  of  the 
error  committed  at  each  range  step,  and  the  cumulative 
range  values  for  the  transmission  loss  file. 

II.  ENHANCED  IMPLEMENTATION  (EIFD) 

Later  in  this  paper,  we  will  vary  the  range  step  so  that 
(A  r)m  =r„  —  rm  _  ,  refers  to  the  mth  range  step  taken  in 
the  calculation.  In  fact,  for  different  problems,  the  behavior 
of  (A r)m  can  actually  vary  in  interesting  ways,  but  often 
(A r)„  can  be  surprising  larger  than  a  wavelength.  Conse¬ 
quently,  significant  computational  advantages  are  obtained. 
Since  we  intend  to  select  the  range  step  size  to  control  a 
measure  of  the  error,  we  first  require  a  method  for  estimating 
the  error  associated  with  a  given  step  size.  One  such  choice  is 

£(rm)  =  (Ar)2m||<3„u;||,,  (5) 

which  indicates  local  error  estimate  and  feedback  to  adjust 
the  range  step  size.  To  estimate  drr  u"  pointwise  by  Eq.  ( 5 )  at 
each  depth  node  n  =  l,...,iV,  we  use  backward  differences  on 
the  two  previous  range  steps;  i.e., 


The  norm  used  in  Eq.  ( 5 )  is  the  discrete  approximation  to 
the  L2  norm,  given  by 


where  the  overbar  denotes  complex  conjugation.  We  exam¬ 
ined  other  norms,  £,  or  L  x  ,  but  found  no  particular  reason 
to  select  one  over  another,  so  that  L2  was  chosen  as  a  matter 
of  convenience.  For  the  PE,  the  measure  given  by  Eq.  (5)  is 
only  proportional  to  the  true  error,  but  is  nonetheless  a  good 
error  indicator.  A  more  accurate  error  estimation  technique, 
such  as  Richardson  extrapolation,  could  also  be  implemen¬ 
ted  at  the  expense  of  additional  computational  overhead. 141 5 
Furthermore,  this  error  indicator  does  not  contain  any  esti¬ 
mate  of  the  error  caused  by  the  discretization  of  the  depth 
into  steps  of  size  Az. 

Our  technique  contrasts  with  the  work  of  Ref.  1 6,  which 
does  not  deal  with  step-changing  algorithms,  but  which  uses 
a  measure  of  the  local  error  to  verify  the  solution  accuracy 
after  the  computation  of  the  entire  solution.  Furthermore, 
we  emphasize  that  our  error  indicator  is  obtained  directly 
from  the  computed  solution,  a  substantially  different  ap¬ 
proach  from  the  oue  used  in  Ref.  3,  which  estimates  trunca¬ 
tion  error  of  the  solution  in  terms  of  certain  partial  deriva¬ 
tives  of  the  refractive  index.  This  latter  method  of  error 
estimation  is  feasible  for  calculations  performed  in  the  wave- 
number  domain. 

We  seek  to  control  the  magnitude  of  the  error  indicator 
Fq  (5)  in  two  ways.  First,  the  error  indicator  E  should  be 
kept  below  a  user-specified  error  tolerance  e  in  order  to 
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maintain  the  solution  accuracy.  Second,  E  should  not  fall 
below  some  percentage,  chosen  as  70,  of  e  so  that  reasonable 
computational  efficiency  is  achieved.  When  these  conditions 
hold,  the  computation  continues  with  the  same  value  of 
( Ar) m  Otherwise,  if  E >  e  or  E  <  0.7e,  the  error  indicator  is 
controlled  by  adjusting  the  step  size  ( Ar )  m .  Our  technique  is 
similar  to  the  step  size  selection  used  by  one-step  codes  for 
ordinary  differential  equations  described  in  Ref.  1 7.  Since  we 
are  using  an  implicit  scheme,  a  system  of  equations,  whose 
matrix  representation  is  tridiagonal,  is  solved  at  each  range 
step.  This  calculation  is  efficient  because  of  the  special  struc¬ 
ture  of  the  matrix.  Since  elements  in  this  system  are  depen¬ 
dent  on  (A r)m,  it  is  desirable  to  avoid  modifications  of 
( Ar)m  at  every  step  and  the  costly  recomputation  of  matrix 
elements.  We  have  found  that  a  proper  trade-off  occurs  with 
the  indicator  inside  a  tolerance  window  of  from  70%  to 
100%  of  the  error  tolerance. 

After  any  solution  step,  the  next  range  step  ( A r)m  .  ,  is 
determined  from  E,  e  and  (A r)m.  For  the  current  step,  we 
know  from  Eq.  ( 5 )  that  E  a  ( Ar )  jj,  and,  to  achieve  accurate 
and  efficient  computation,  we  want  the  next  range  step  size 
( A r)m+  ,  a  e l/:.  The  following  proportionality  results: 

e/£(rm)oc(Ar)i+1/(Ar)Jm.  (8) 

Consequently,  if  a  range  step  adjustment  is  judged  necessary, 
the  new  step  size  is  determined  by 

(Ar)m+,  =  £  [e/£(rm )  ]  l/2(Ar)m,  (9) 

where  £  is  a  constant  of  proportionality. 

The  parameter  £  provides  a  useful  degree  of  freedom  in 
our  method.  If  the  current  step’s  error  estimate  is  over  the 
window  (£>  e),  the  adjustment  is  made  with  £  =  fOJ  in 
order  to  project  the  error  estimate  for  the  next  step  toward 
the  bottom  of  the  window.  If  the  error  estimate  has  de¬ 
creased  below  the  window  (£<0.7e),  the  choice  £  =  1  ad¬ 
justs  the  next  step  size  to  an  error  estimation  near  the  top  of 
the  window.  These  selections  of  £  provide  for  maximum 
possible  use  of  the  window  before  the  next  range  step  adjust¬ 
ment  needs  to  be  made.  We  emphasize  the  novelty  of  our 
approach  by  pointing  out  that  selection  of  parameters  analo¬ 
gous  to£  are  usually  ad  hoc . 17  In  contrast,  we  have  provided 
herein  a  definite  and  consistent  way  to  select  its  value  based 
on  the  performanci  f  v  error  indicator. 

Another  type  of  ?d  rrment  to  ( Ar)„  +  ,  may  also  be 
required  to  ensure  tha  algorithm  performs  adequately  at 
vertical  interfaces,  which  may  occur  at  certain  ranges  to 
model  rapid  horizontal  changes  in  sound-speed  profile,  bot¬ 
tom  structure,  or  channel  geometry.  As  a  result  of  the  basic 
IFD  implementation,  it  is  necessary  that  calculation  always 
begin  exactly  at  a  range  r„,  of  a  vertical  interface.  With  uni¬ 
form  step  sizes,  this  is  usually  easy  to  arrange,  but  with  vari¬ 
able  range  steps,  particularly  when  (Ar)m  is  very  large,  it  is 
possible  for  (rv,  -  rm  )  « ( Ar)„,.  In  this  case.  (A r)m„, 
=  rvl  —  rm  and  the  range  step  is  forced  to  be  much  smaller 
than  required  by  the  error  indicator.  Once  the  computation 
proceeds  past  the  vertical  interface,  many  more  step  adjust¬ 
ments  arise  from  this  small  step,  introducing  undesirable 
computational  overhead.  In  order  to  provide  a  smooth  tran¬ 
sition  of  the  range  step  size  at  such  an  interface,  our  algo- 
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rithm  checks  two  step  sizes  ahead  and,  when  a  change  is 
necessary,  makes  two  modest  step  size  changes  instead  of 
one  more  drastic  change.  In  particular, 

if  rm  +  2(A r)m+  ,  >rv, 

then  =  $(rv,  -  rm).  (10) 

Occasionally,  the  step  size  from  Eq.  ( 10)  is  too  large  and  the 
indicator  leaves  the  window.  However,  within  two  steps,  the 
error  is  typically  brought  back  under  control. 

EIFD  includes  another  enhancement  appropriate  for 
acoustic  intensity  calculations.  Any  algorithm  that  numeri¬ 
cally  solves  a  PE  does  not  compute  the  wave  field  but,  rather, 
the  slowly  modulated  envelope  u(rj).  This  quantity  is  then 
used  to  calculate  transmission  loss  or  relative  intensity. 
Thus,  to  underwater  acousticians,  a  meaningful  error  toler¬ 
ance  would  be  expressed  in  decibels,  a  form  of  relative  error, 
while  a  more  natural  way  to  control  error  is  in  terms  of  the 
absolute  units  in  which  the  envelope  function  is  computed. 
Clearly,  near  regions  where  the  solution  norm  is  zero  or 
nearly  so,  the  decibel  error  could  be  very  large  and,  in  fact, 
may  be  impossible  to  control.  Thus  any  adaptive  method 
should  be  constructed  to  locate  deep  fades  and  indicate  their 
size  to  some  extent,  but  not  to  accurately  predict  how  many 
decibels  are  lost  at  the  fade.  We  now  describe  one  way  to 
relate  the  absolute  error  e  to  a  relative  error  q,  expressed  in 
decibels. 

Let  u  be  the  true  solution  to  Eq.  ( 3 )  and  u  the  computed 
solution  at  some  mesh  point,  with  £  the  computed  value  of 
the  error  indicator.  We  seek  an  estimate  for  the  error  q  in 
decibels.  Note  that 

||fi|  -  \u\\<kE,  (11) 

where  k  is  a  positive  constant.  Equation  (11)  is  one  way  to 
state  that  the  true  error  is  proportional  to  the  error  indicator. 
This  expression  can  be  reformulated  as 

20  logl0(  1  -  /cJE" /j w; ) <20  Iog,„|w|  -  20  log10|ti| 

<201oglo(l  +kE/\u\).  (12) 

From  Eq.  ( 12),  we  get  bounds  on  the  decibel  error  q,  specifi¬ 
cally 

q  =  |20  logloju|  -20  loglo|u|| 

<20|loglo(  1  —  kE /\u\)\.  (13) 

In  general,  q  is  small  since  the  error  tolerance  e  is  typically 
selected  to  force  £<  |  u  \ .  However,  near  a  null  or  deep  fade,  q 
could  become  large  since  E  may  equal  or  substantially  ex¬ 
ceed  { u | .  This  is  one  reason  why  the  transmission  loss  (or 
relative  intensity),  which  is  measured  in  decibels,  when  cal¬ 
culated  adaptively,  may  exhibit  substantial  error  in  the  vi¬ 
cinity  of  deep  fades  and  other  locations  where  the  solution 
magnitude  and  error  tolerance  are  of  comparable  size. 

We  can  use  Eq.  (13)  to  generate  a  choice  of  error  toler¬ 
ance  for  a  prescribed  level  of  decibel  error  q,  which  will  be 
valid  away  from  deep  fades.  It  follows  from  Eq.  (13)  and 
£>  0.7e  that 

e<(10|u|/7*)(l  -  10-’/2°).  (14) 

Air  optimal  numerical  value  for  k  i»  uncertain,  but  from  our 
numerous  computations,  it  appeared  that  kz:  0.1  in  many 
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cases.  An  uncertainty  of  q  =  1  dB  seems  to  be  a  reasonable 
constraint  on  the  size  of  the  decibel  error  away  from  deep 
fades.  Given  these  choices,  a  convenient  selection  for  a  value 
of  the  error  tolerance,  using  Eq.  ( 14),  is  e  =  10“ 3.  Other 
choices  of  e  are  discussed  in  Sec.  III. 

Before  discussing  specific  numerical  examples,  we  now 
discuss  several  additional  features  that  have  been  included  in 
EIFD.  Often  the  start-up  field  results  in  many  modes  propa¬ 
gating  until  bottom  attenuation  strips  the  higher-order 
modes  so  that  near  the  source,  the  solution  looks  “noisy.” 
Within  this  region,  the  adaptive  range  step  changer  could 
make  the  range  step  very  small.  Since  the  nearfield  is  rarely 
of  interest  in  calculations  with  EIFD,  a  feature  has  been 
added  to  the  input  runstream  that  allows  the  adaptive  calcu¬ 
lation  to  be  switched  on  after  a  certain  range,  called  the  cut- 
in  range  ra,  is  achieved.  For  example,  in  a  particular  ocean 
channel,  it  may  be  known  that,  using  the  usual  Gaussian 
start-up  field,  the  effects  of  bottom  attenuation  are  largely 
absent  until  a  range  of  4  km.  Consequently,  ra  is  set  to  this 
value. 

It  is  sometimes  possible  that  the  solution  to  the  PE  is 
exponentially  decaying.  In  this  circumstance,  the  solution 
norm  will  eventually  become  smaller  than  the  error  toler¬ 
ance,  so  that  the  error  measure  becomes  useless.  EIFD  con¬ 
tains  a  feature  that  checks  the  norm  of  the  solution  vector  at 
each  step.  If  the  norm  ||u||:  is  less  than  ten  times  the  error 
tolerance  e,  then  the  step  size  is  no  longer  changed. 

In  order  to  monitor  the  progress  of  the  adaptive  calcula¬ 
tion,  several  additional  output  files  are  required.  Specifical¬ 
ly,  both  the  error  indicator  and  the  step  size  at  each  step  are 
listed  to  output  files.  Furthermore,  the  value  of  each  range  at 
which  the  relative  intensity  is  written  must  also  be  listed,  so 
that  intensity  curves  can  be  properly  plotted. 

III.  NUMERICAL  RESULTS 

In  order  to  gauge  the  efficacy  of  our  algorithm,  four 
examples  will  be  discussed,  with  each  example  presenting  a 
different  and  acoustically  important  environment.  We  have 
thoroughly  and  carefully  examined  the  performance  of  our 
step  changing  algorithm,  with  attention  focused  both  on  so¬ 
lution  accuracy  as  well  as  on  execution  time  improvements. 
Although  other,  nonfinite-difference  implementations  have 
utilized  range  step  changers,  an  examination,  like  ours,  of 
algorithm  performance  is  evidently  novel. 

All  calculations  were  performed  on  a  Prime  850  mini¬ 
computer.  Efficiences  for  each  example  are  estimated  using  a 
performance  measure  that  is  the  run-time  ratio  of  the  nona- 
daptive  (.V  mode)  run  time  to  the  adaptive  (A  mode)  run 
time.  This  ratio,  denoted  by  the  symbol  £.  suggests  the  gen¬ 
eral  increase  in  computational  efficiency  of  the  A  mode  when 
£>  1,  but  is  naturally  dependent  to  some  extent  on  the  hard¬ 
ware  used.  Nonetheless,  we  believe  that  £is  both  a  straight¬ 
forward  way  to  measure  the  performance  of  our  algorithm 
and  a  reasonable  indicator  of  speed  ratios  on  any  system. 

In  each  of  our  examples,  the  calculated  quantity  is  rela¬ 
tive  intensity  I,  given  by 

I  =  20  log  10i  (15) 

where p(r,z)  is  given  by  Eq.  ( 2 )  and  prtf  is  a  reference  pres- 
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sure  measured  at  r  =  1  m  from  the  source.  Furthermore,  in 
each  example,  a  Gaussian  starting  field  is  used  to  provide  an 
estimate  of  the  sound  field  close  to  the  source.  This  method  is 
widely  used  because  it  is  simple  to  calculate,  and,  in  fact,  is 
an  integral  part  of  the  standard  implementation  of  IFD. 
Nonetheless,  other  starting  methods  are  available  and  could 
be  used  here  as  well.  However,  our  results  are  not  expected  to 
depend  significantly  on  the  choice  of  a  Gaussian  starting 
field. 

A.  Isospeed  horizontal  channel 

In  this  example,  a  cw  source  S  of  frequency  200  Hz  is 
placed  at  depth  25  m.  The  water  is  100  m  deep,  overlying  a 
fluid  bottom  that  extends  an  additional  100  m.  Finally,  an 
artificial  absorbing  layer  of  depth  50  m  is  appended.  Sound 
speed  and  density  discontinuities  appear  at  the  interface,  and 
the  lower  layer  also  has  volume  attenuation.  The  channel 
geometry,  together  with  values  of  important  acoustical  pa¬ 
rameters,  is  shown  in  Fig.  1.  This  sound  channel  is  identical 
to  one  used  in  Ref.  10. 

For  this  channel,  the  adaptive  calculation  was  accom¬ 
plished  by  selecting  an  error  tolerance  of  e  =  10  - '  and  cut- 
in  range  of  ra  =  3  km.  The  behavior  of  the  error  indicator  E 
versus  range  r  is  shown  in  Fig.  2.  The  step  size  begins  at 
Ar  =  5  m,  and  remains  constant  until  the  cut-in  range  is 
achieved.  The  error  indicator  decreases  steadily  as  the 
“noisy"  start-up  field  is  stripped  by  the  lossy  bottom.  At  that 
point,  E  is  below  the  tolerance  window,  indicated  by  the  pair 
of  horizontal  dashed  lines.  Once  the  step  size  changer  is  acti¬ 
vated,  A r  is  increased  and  E  is  brought  inside  the  window. 
The  error  continues  to  decrease,  and  A r  is  enlarged  each  time 
E  falls  below  the  lower  bound  of  the  window.  At  about  range 
r  =  36  km,  the  step  size  is  approximately  220  m,  or  about  26 
wavelengths.  That  the  step  size  could  achieve  such  a  large 
size  and  still  retain  accuracy  of  calculation  may  be  surpris¬ 
ing. 

The  accuracy  of  the  adaptive  calculation  is  clearly  dem¬ 
onstrated  in  Fig.  3.  The  solid  curve  represents  relative  inten¬ 
sity  /  versus  range  r  for  a  receiver  depth  of  25  m,  using  N 
mode  with  Ar  =  5  m.  The  dashed  curve  is  the  same  quantity 
but  computed  in  A  mode  as  described  above.  The  overall 
shape  of  the  intensity  curve,  including  the  location  of  local 
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FIG.  I.  Isospeed  horizontal  channel. 
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FIG.  2  Error  indicator  E  versus  ranger  for  channel  shown  in  Fig.  1  :/=  200 
Hz;  h,  =  25  m;  r,  =  5  km;  and  e  =  10" ' 


maxima  and  minima  and  height  of  the  peaks,  is  predicted 
well  by  the  adaptive  calculation.  The  extent  of  deep  fades  can 
be  over-  or  underestimated  by  several  decibels  or  more.  Al¬ 
though  this  is  expected  from  our  earlier  discussion  in  Sec.  II, 
it  is  of  little  consequence  since  intensity  levels  at  deep  fades 
are  usually  of  little  practical  interest.  For  the  example,  the 
efficiency  factor  is  F  =  7.3,  meaning  that  the  adaptive  calcu¬ 
lation  ran  over  seven  times  faster  than  the  nonadaptive  one. 

B.  Converging  channel 

In  this  example,  the  sound  channel  is  range  dependent. 
As  depicted  in  Fig.  4,  an  isospeed  water  column  with 
c0  =  1500  m  s~ 1  begins  at  depth  350  m  and  lies  over  an 
interface,  with  a  fluid  layer  below  that  possesses  volume  at¬ 
tenuation.  The  interface  is  horizontal  for  a  range  of  10  km. 
Then,  the  interface  slopes  upward  at  an  angle  of  8.5  deg  until 
range  12  km,  where  the  depth  is  50  m.  At  that  point,  the 
interface  is  flat  and  remains  so  to  a  distance  of  40  km.  The 
lower  layer  always  extends  to  depth  750  m,  and  an  artificial 


FIG.  3  Relative  intensity  /  versus  range  rfor  parameters  of  Fig.  2;  h,  =  25 
m. 
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FIG  4  Isospeed  sloping  channel 

absorbing  layer  extends  an  additional  250  m.  This  channel 
was  described  in  Ref.  13.  The  source  frequency  is  25  Hz  and 
is  located  at  depth  25  m. 

The  error  indicator  as  a  function  of  range  is  depicted  in 
Fig.  5.  The  step  size  was  held  fixed  at  Ar  =  10  m  until  the 
cut-in  range  ru  =  4  km.  Note  that  the  error  tolerance  used  in 
this  example  is  e  =  1.5 x  10~4.  This  value  was  chosen  to 
force  the  tolerance  and  the  error  indicator  to  have  approxi¬ 
mately  similar  magnitudes  at  the  cut-in  range,  preventing 
rapid  and  severe  step  size  adjustments.  The  step  size  contin¬ 
ues  to  increase  until  Ar  =  25  m  at  the  foot  of  the  slope  where 
r  =  1 0  km.  The  sudden  rapid  change  in  E  seen  at  that  point  is 
caused  by  step  truncation  at  the  vertical  interface  at  the  start 
of  the  slope  as  was  discussed  in  Sec.  II.  As  the  sound  propa¬ 
gates  upslope,  the  error  starts  to  increase  steadily  and  is  con¬ 
trolled  by  reducing  step  size,  until  Ar  =  9  m  at  the  top  of  the 
slope.  Again,  a  step  truncation  occurs  at  r  =  12  km  for  the 
reason  just  mentioned.  After  E  is  under  control,  the  error 
decreases  in  a  somewhat  regular  manner,  permitting  step 
size  to  be  increased  until  Ar  =  96  m  at  range  near  r  =  24  km. 
This  step  size  corresponds  to  1.6  wavelengths.  Beyond  24 
km,  E  decreases  rapidly  and  no  more  step  adjustments  are 
made.  The  reason  is  that  the  solution  norm  is  too  small  and 
the  algorithm  has  stopped  changing  step  size.  For  this  par¬ 
ticular  sound  channel,  the  step  size  changer  makes  nearly 
full  use  of  the  tolerance  window.  We  note  that  this  error 


FIG  5  Error  indicator  E  versus  range  r  for  the  isospeed  sloping  channel 
shown  in  Fig  4.  with  source  /  at  the  left  and  receiver  .ft  at  the  right./  =23 
Hi.  h  =  25  m:  r„  =  4  km.  and  (  =  I  5  x  I0"J. 


FIG  6.  Relative  intensity /versus  range  r  for  parameters  of  Fig.  5:  A.  =  25 
m. 

behavior,  while  desirable  in  every  case,  is  not  fully  achieved 
in  all  sound  channels  (see  Fig.  2). 

A  comparison  between  N  mode  and  A  mode  intensities 
is  shown  in  Fig.  6.  As  in  our  first  case,  the  receiver  depth  is  25 
m.  The  solid  curve  depicts  relative  intensity  versus  range  for 
the  N  mode  solution,  while  the  dashed  curve  depicts  1  for  the 
A  mode  calculation.  As  can  be  seen,  acoustic  mode  cutoff 
occurs  near  the  top  of  the  slope.  Thus  the  sound  field  decays 
exponentially  beyond  this  range,  accounting  for  why  the  step 
changer  stopped.  Some  error  is  noticeable  between  the  two 
intensity  curves,  but  the  agreement  is  extremely  good.  Fur¬ 
thermore,  the  efficiency  factor  F  —  2.3.  While  not  as  dra¬ 
matic  as  the  isospeed  problem,  this  result  represents  a  sub¬ 
stantial  improvement  in  computation  time. 

C.  Diverging  channel 

This  channel  is  physically  identical  to  the  one  just  de¬ 
scribed.  However,  the  source  J  and  receiver  ^  have  been 
interchanged,  so  that  the  sound  propagates  to  10  km  in  very 
shallow  water,  travels  down  a  steep  slope,  then  propagates 
out  to  40  km  in  water  350  m  deep.  Pertinent  environmental 


FIG  7  Error  indicator  E  versus  range  r  for  channel  shown  in  Fig  4.  but 
with  source  /  and  receiver  .w  interchanged,  parameters  as  in  Fig.  6 
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parameters  and  source/ receiver  depths  remain  the  same. 

The  behavior  of  the  error  indicator  is  shown  in  Fig.  7. 
Values  for  the  error  tolerance  and  cut-in  range  are  the  same 
as  those  used  in  Fig.  5.  Beyond  ra,  step  size  is  rapidly  in¬ 
creased  until  Ar  =  42  m  at  the  top  of  the  slope.  After  the 
occurrence  of  step  truncation  at  that  point,  E  decreases  each 
time  A r  is  changed.  Consequently,  A r  is  increased  as  the 
sound  propagates  downslope,  which  is  the  opposite  behavior 
from  what  occurred  when  propagation  was  upslope.  At 
r  =  12  km,  there  is  another  step  truncation  and  Ar  changes 
increases  times.  Near  14  km,  a  marked  decrease  occurs  in  the 
rate  of  decrease  of  E  with  increasing  range.  In  fact,  the  final 
change  in  step  size  occurs  at  r  —  24  km,  where  A r  =  168  m, 
corresponding  to  2.8  wavelengths.  This  step  size  is  substan¬ 
tially  larger  than  the  terminal  step  used  in  Fig.  5. 

Computed  intensities  are  depicted  in  Fig.  8.  The  solid 
curve  represents  relative  intensity  /  without  adaption,  while 
the  dashed  curve  shows  intensity  calculated  with  the  error 
tolerance  and  cut-in  range  of  Fig.  7.  Agreement  is  excellent 
between  both  calculations  up  through  and  well  past  the  slop¬ 
ing  region.  Beyond  about  r  —  20  km,  there  are  some  observ¬ 
able  phase  shifts  in  both  peaks  and  fades.  Nonetheless,  the 
qualitative  picture  of  fades  and  peaks  intensity  in  the  chan¬ 
nel.  along  with  levels  at  most  of  the  peaks,  is  represented  very 
well.  Recall  that  in  the  previous  example,  adaptation  was 
terminated  at  r  =  24  km,  where  the  relative  intensity  was 
about  —  100  dB.  In  Fig.  8,  the  solution  norm  beyond  the 
slope  ts  just  barely  above  the  termination  criterion.  Conse¬ 
quently,  inaccuracies  inevitably  creep  into  the  computation. 
Naturally,  more  accuracy  can  be  attained  with  stricter  val¬ 
ues  of  the  error  tolerance  e.  Finally,  the  efficiency  factor  is 
F  —  4.8  for  the  two  runs  of  Fig.  8.  This  figure  is  substantially 
higher  than  that  for  Fig.  6  for  the  same  physical  channel,  and 
indicates  a  nearly  fivefold  decrease  in  run  time  for  the  A 
mode  calculation. 

D.  Deep-ocean  channel 

In  this  example,  a  deep-ocean  channel  with  moderately 
refracting  sound-speed  profile  was  used.  This  example  is 
similar  to  one  used  in  Ref.  18.  The  profile  contains  a  front. 


RAN6E  (KM) 


FIG.  9.  Level  curves  of  sound  speed  (ms*1)  versus  range  and  depth  ( km ) 
for  a  deep-water  sound  channel. 


located  at  range  50  km,  as  suggested  in  Fig.  9.  The  sound 
speed  is  given  by 

c(rj)  =c„[l  +^{r)(e-’>  +  v-  D],  (16) 

where  ca,  g(r),  and  rj(r,z)  are  given  in  Table  I.  This  sound- 
speed  profile  is  a  modification  of  Munk’s  canonical  deep¬ 
water  profile.  The  channel  is  5  km  deep  and  an  artificial 
absorbing  bottom  of  depth  1  km  is  used.  The  source  is  at 
depth  100  m,  and  the  source  frequency  is  100  Hz. 

The  behavior  of  the  error  indicator  is  shown  in  Fig.  10. 
For  this  calculation,  the  error  tolerance  is  e  =  10~3  and  cut- 
in  occurs  at  ra  =40  km.  The  cut-in  distance  is  large  because 
the  error  indicator  E  decreases  much  more  slowly  than  in 
any  of  the  previous  three  examples,  owing  to  the  existence  of 
a  sound  channel  and  a  lack  of  volume  attenuation  in  the 
bottom  (the  artificial  absorbing  layer  prevents  reflections 
from  the  bottom  boundary  but  does  not  induce  exponential 
decay  in  the  solution).  The  initial  step  size  is  A r  =  10  m, 
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FIG  8  Relative  intensity  /  versus  range  rfor  parameters  in  Fig.  6. 


TABLE  I.  Components  and  parameters  for  >he  range-dependent  sound- 
speed  profile  of  Eq.  (16). 


Parameter  Value  or  description 

c„  sound  speed  at  channel  axis,  1450  ms"1 

f(r)  B(r)g/c . 

B{r)  thickness  of  thermocline  front, 

B ,  +  [(flj  -  B|)/2)tanhl(r—  r,\/L\ 

S  1.2  km 

B-_  1 .0  km 

rf  range  to  front  center,  60  km 

L  front  width,  20  km 

g  0.017 

rj(r.z)  2{[z  -  Z,(r)]/B(r)} 

z„  ( r)  depth  of  channel  axis, 

fu  —  I(ij  —  B,)/2)tanh((r-  rr)/L] 

initial  depth  of  sound  channel  axis  at  r  =  0  km,  0.8  km. 
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FIG.  10  Error  indicator  £  versus  range  r  for  the  environment  of  Fig.  9. 
/=  100  Hz:  h .  =  100  m;  r„  =  40  km:  and  f  =  10 ' '. 


which  was  used  in  Ref.  18.  At  cut-in,  the  error  indicator 
happens  to  be  just  outside  the  window  and  the  step  size  drops 
to  Ar  =  7  m.  At  about  r  =  42  km,  the  error  has  fallen  enough 
so  that  the  first  step  size  increase  occurs.  Then,  several  step 
increases  in  a  row  occur.  Note  that  the  error  indicator  does 
not  increase  significantly  in  this  region.  Then,  at  nearly  50 
km,  E  becomes  virtually  constant  for  almost  1 5  km.  Finally, 
it  leaves  the  window,  the  step  size  is  increased,  and  this  time 
the  error  indicator  is  near  the  center  of  the  window.  At 
r  =  105  km,  the  step  size  is  increased  again  to  A r  =  49  km, 
and  the  error  is  seen  to  remain  virtually  flat  out  to  the  maxi¬ 
mum  range  examined  in  this  example.  It  would  appear  that, 
in  this  last  region,  the  algorithm  has  found  a  nearly  optimal 
step  size  (about  three  wavelengths),  in  the  sense  that  the 
error  is  virtually  constant  and  the  error  is  near  the  center  of 
the  tolerance  window. 

Relative  intensity  curves  for  this  sound  channel  are 
shown  in  Fig.  1 1 .  The  solid  curve  depicts  /  versus  range  r  for 
a  receiver  depth  of  300  m  computed  in  N  mode,  while  the 


FIG.  II,  Relative  intensity  /  versus  range  r  for  the  deep-ocean  channel 
shown  in  Fig.  9;  A,  «  300  m;  other  parameters  as  in  Fig.  10. 


dashed  curve  represents  the  intensity  for  the  A  mode  calcula¬ 
tion.  The  two  curves  are  in  excellent  agreement  throughout 
the  first  convergence  zone,  between  about  r  =  45  km  to 
r  =  70  km,  and  out  through  part  of  the  second  convergence 
zone.  Beyond  about  r  =  1 1 5  km,  there  appears  to  be  some 
shifting  of  phase  pattern.  The  A  mode  curve  reflects  the 
qualitative  behavior  of  intensity  all  the  way  to  r  =  140  km. 
For  this  example,  the  efficiency  factor  F=  1.8.  One  reason 
for  the  lower  value  of  F  is  the  size  of  the  cut-in  range  ra 
While  smaller  than  F  values  attained  for  the  previous  three 
cases,  this  still  represents  a  substantial  savings  of  processor 
time,  especially  for  this  computationally  intensive  example. 

IV.  SUMMARY 

The  parabolic  approximation  to  the  reduced  wave  equa¬ 
tion  has  established  itself  as  a  formidable  propagation  model 
within  the  underwater  acoustics  community.  Because  it  in¬ 
cludes  range-dependent  environments,  is  valid  at  low  fre¬ 
quencies,  and  can  be  numerically  solved  with  efficient  algo¬ 
rithms,  it  has  become  the  model  of  choice  for  many 
applications.  After  very  briefly  reviewing  the  origin  of  the 
parabolic  approximation  and  related  parabolic  equations, 
we  summarize  pertinent  details,  such  as  notation  and  error 
characteristics  of  the  widely  employed  IFD  implementation. 

We  introduce  one  way  to  estimate  the  error  from  the 
IFD  algorithm.  This  error  only  accounts  for  discretization 
in  range.  The  error  indicator  is  calculated  with  an  appropri¬ 
ate  norm  at  each  step.  For  a  specified  error  tolerance,  a  toler 
ance  window  is  formed  by  requiring  that  the  error  indicator 
remain  below  the  error  tolerance  but  above  70%  of  that  val¬ 
ue.  If  the  error  indicator  leaves  the  window,  the  step  size  is 
either  increased  or  decreased.  The  magnitude  of  the  change 
is  computed  from  a  relationship  between  step  sizes  and  er¬ 
rors.  An  additional  feature  is  that  the  activation  of  the  step 
changer  can  be  postponed  by  specifying  an  additional  pa¬ 
rameter  called  the  cut-in  range.  Below  this  range,  no  range 
step  changes  are  made,  which  may  be  desirable  as  the  rapid 
oscillations  sometimes  present  in  the  start-up  field  are 
stripped  away.  Also,  should  the  ratio  of  the  norm  of  the 
solution  and  the  error  tolerance  drop  below  a  specified  value, 
adaption  is  terminated.  In  addition  to  the  standard  input  and 
output  required  by  IFD,  our  enhancements  (collectiyely 
known  as  EIFD)  require  additional  input,  namely  the  error 
tolerance  and  cut-in  range,  and  generate  several  additional 
output  files  useful  for  monitoring  the  error  behavior  and 
graphically  interpreting  the  transmission  loss  or  intensity. 
Because  the  algorithm  controls  absolute  error,  but  not  rela¬ 
tive  error,  the  decibel  error  in  the  vicinity  of  some  fades  may 
occasionally  be  large.  Nonetheless,  our  adaptive  enhance¬ 
ments  locate  these  fades  and  indicate  their  approximate 
sizes. 

Numerical  examples  illustrate  the  performance  of  our 
algorithm  in  a  variety  of  propagation  channels.  In  a  shallow 
isospeed  channel,  the  adaptive  algorithm  produced  nearly  a 
sevenfold  improvement  (decrease)  in  run  time,  due  partially 
to  its  selection  of  large  step  sizes.  Two  cases  of  propagation  in 
an  isospeed  channel  with  a  sloping  bottom  were  examined 
also.  In  the  first  instance,  the  source  was  in  moderately  deep 
water  and  the  receiver  was  in  very  shallow  water,  while  in  the 
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second,  the  source  and  receiver  were  interchanged.  Step 
sizes  were  found  to  exhibit  interesting  changes  for  propaga¬ 
tion  either  up  or  down  the  slope.  In  the  former  case,  mode 
cutoff  occurred  at  the  top  of  the  slope,  eventually  leading  to 
the  requirement  that  adaption  terminate  because  so  ution 
norm  and  error  tolerance  approached  comparable  magni¬ 
tudes.  In  the  latter,  multiple  modes  were  not  excited  until 
well  down  the  slope.  It  is  interesting  to  note  that  the  adaptive 
algorithm  ran  faster  for  the  downslope  example.  Finally, 
deep-ocean  propagation  through  a  front  was  studied.  Effi¬ 
ciency  according  to  our  measure  increased  by  less  than  a 
factor  of  2,  though  this  represents  an  enormous  savings  in 
computational  time. 

We  have  described  enhancements  to  the  IFD  implemen¬ 
tation  of  the  parabolic  approximation  that  can  yield  signifi¬ 
cant  run-time  improvements.  An  additional  approach, 
which  may  also  offer  promise,  is  to  adaptively  add  (delete) 
points  to  ( from )  the  depth  mesh  grid.  An  optimal  method  of 
simultaneously  selecting  both  range-  and  depth-mesh  incre¬ 
ments  to  minimize  solution  error  is  being  investigated.  These 
on-going  developments  suggest  additional  research  direc¬ 
tions  that  can  significantly  extend  the  family  of  underwater 
acoustic  propagation  problems  that  can  be  solved  on  current 
research  computers. 
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f inite-dif f erence (IFD)  implementation  of  the  parabolic  approximation. 

An  error  indicator  is  computed  at  each  range  step,  and  its  value  is 
compared  to  an  error  tolerance  window  that  is  readily  specified  by  the 
user.  If  the  error  indicator  falls  outside  this  window,  a  new  range 
step  size  is  computed  and  used  until  the  error  indicator  again  leaves 
the  tolerance  window.  Furthermore,  for  a  given  tolerance,  this  algorithm 
generates  a  range  step  size  that  is  optimal  in  a  specified  sense  and 
that  often  leads  to  large  decreases  in  run  time.  Additional  related 
modifications  to  the  IFD  implementations  are  discussed.  Several  examples 
are  presented  that  illustrate  the  efficacy  of  the  enhanced  algorithm. 
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